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Abstract 


Several boundary conditions that allow subsonic and supersonic flow into and out of the 
computational domain are discussed. These boundary conditions are demonstrated in 
the FUN3D computational fluid dynamics (CFD) code which solves the three-dimensional 
Navier-Stokes equations on unstructured computational meshes. The boundary conditions 
are enforced through determination of the flux contribution at the boundary to the solution 
residual. The boundary conditions are implemented in an implicit form where the Jacobian 
contribution of the boundary condition is included and is exact. All of the flows are gov- 
erned by the calorically perfect gas thermodynamic equations. Three problems are used to 
assess these boundary conditions. Solution residual convergence to machine zero precision 
occurred for all cases. The converged solution boundary state is compared with the re- 
quested boundary state for several levels of mesh densities. The boundary values converged 
to the requested boundary condition with approximately second-order accuracy for all of 
the cases. 

1 Introduction 

T he solution of a partial differential equation or a system of partial differential equations 
requires a statement of the dependent variables along the boundaries of the solution 
space. When physically relevant solutions of the equations are sought, the imposition of the 
boundary conditions should solve the correct numerical problem without adversely affecting 
solution stability. 

The choice between an explicit or implicit implementation of the boundary conditions 
in the solution algorithm is dependent on the numerical method that is used for solving the 
governing equations. Explicit treatments of boundary conditions has occupied much of the 
literature. The compilation of papers in the proceedings Numerical Boundary Condition 
Procedures [1] contains extensive discussions of the state-of-the-art theories at the time of 
printing. In particular, Moretti [2], Blottner [3], Pulliam [4], and Thompkins and Bush [5] 
discuss inflow and outflow boundary conditions for both Navier-Stokes and Euler equations. 

Much of that work centered around developing a system of equations, often in terms of 
characteristic variables, to describe the desired physical state in a consistent form at the 
boundary. Information is communicated through the boundary interface via wave propaga- 
tion and convection. The flow of information across the boundary determines the physical 
conditions to impose. Thus, the boundary conditions must be formulated to keep the so- 
lution realizable at the boundary. The interior solution will be consistent with the speci- 
fied physical state at all of the boundaries once an iteratively converged solution has been 
achieved. 

Figure 1 is a schematic of a boundary where the physical state is on the right side and 
the solution space is on the left. Also shown are the upstream and downstream propagating 
Riemann invariants, R~ and R + . The convention in FUN3D [ 6 ] is that the normal vector at 
a boundary, n = n/|n|, is outwardly directed from the interior. The velocity perpendicular 
to the boundary, U_l, is defined by the scalar product of the local velocity vector with the 
boundary normal, U_l = U • n, where U = ui + vj + wk. 

Designated by the sign and magnitude of U_l, the flow conditions at the boundary are 
listed in Table 1. The Euler equations have five eigenvalues — three that are associated with 
convective waves, A 2 - 4 , and two that are associated with acoustic waves, Ai^. A positive 
eigenvalue corresponds to a wave that is entering the domain and that conveys physical 
information specified from the outside (i.e., the boundary condition). A negative eigenvalue 
is a wave that is leaving the solution domain. As discussed in Kim et al. [7], various 
conditions for a boundary are summarized in Table 1. For example, for the specification 
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Figure 1. Characteristics at a boundary. 


of a subsonic inflow boundary with the use of primitive variables, density and velocity or 
density and pressure are combinations that can be used to completely define the condition. 
In practice, many methods are used to enforce a particular condition at a boundary interface. 
In this work, the boundary conditions are enforced weakly, that is, by calculating and 
applying the flux contribution of the boundary condition to the solution residual rather 
than setting the physical value of the boundary condition. 

The details of each boundary condition are discussed in a generic framework in section 2; 
and a discussion of the discrete, computational method used for evaluation of the boundary 
conditions in FUN3D follows in section 3. Three test geometries are described in section 4. 
Computations that demonstrate the use of the boundary conditions are given in section 5. 

Table 1. Flow Conditions at a Boundary 


Case Condition 


Eigenvalue 






Ai 

A 2 

A 3 

A 4 

A5 


1 

Subsonic inflow 

<0 

<c 

<0 

>0 

>0 

>0 

>0 

Pi, T t, a , P* 

2 

Subsonic outflow 

>0 

<c 

<0 

<0 

<0 

<0 

>0 

P 

3 

Supersonic inflow 

<0 

>c 

>0 

>0 

>0 

>0 

>0 

All 

4 

Supersonic outflow 

>0 

>c 

<0 

<0 

<0 

<0 

<0 

None 


* Alternatively p and U or p and p, but not U and p. 


2 Boundary Conditions 


To evaluate the fluxes at the boundary, information from the interior of the domain is 
combined with the physical constraints of the problem. As shown in Figure 2, the flux at 
the boundary depends on a combination of the left and right states (q L and q R ) and the 
direction, n, and area, A, of the boundary (fF = f(n, A, q L , q R )). The left state is equal to 
the interior state, q ( , and the right state is a function of the interior state, q b the free-stream, 
q oal and/or user-specified parameters, 3, depending on the boundary condition. 

The linearization of the flux function for fixed grids, is written as the change in the flux 
function, F ', with respect to the state vector, q , on the left side of the boundary: 


OF dT dF dg R 

dq L dq L + dq R dq L 
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Right state 


Q l = Qi 


<?R = 


Figure 2. Flux and boundary states. 


The first two derivatives, d^F/dq L and dT/d q R , are the Jacobians for the numerical flux 
function used in the spatial differencing. FUN3D provides multiple options for the numerical 
flux function. The results in this work use Roe’s approximate flux scheme and the Jacobian 
terms of Roe’s scheme were handcoded term-by-term. The third term, dq R /dq L , is deter- 
mined exactly by using automatic differentiation via operator overloading of the equations 
in the boundary condition functions in Griewank and Corliss [8,9]. This differentiation 
technique, which is based on the chain rule, permits the partial derivatives of the boundary 
condition variables with respect to the interior variables to be calculated and propagated 
through each of the boundary condition function calls and results in the Jacobian term of 
the boundary conditions. 

Depending upon the specific application, the four inflow/outflow conditions listed in 
Table 1 ( that is, subsonic inflow/outflow and supersonic inflow/outflow ) can be formed in 
different ways. Several methods for determining the boundary state, q R , are discussed in 
sections 2.1 through 2.9. For each boundary condition, the format of the introductory table 
is as follows: the type of condition, the specified conditions, variables that are extrapolated 
from the interior, and variables that are updated as a consequence of applying the boundary 
condition. Words expressed with the fixed-width font indicate code variables that occur 
within FUN3D. 

Variables with dimensions are denoted by an overhead ~ symbol. Non-dimensional vari- 
ables are denoted without any symbol overhead. Freestream values are denoted by the 
subscript oo . The equations relating the dimensional to the non-dimensional flow-field 
variables are as follows: 


The equations in the subsequent sections are written in their non-dimensional form. 

2.1 Far-Field Boundary Condition (f arf ield_roe) 


The far-held boundary condition specifies the free-stream conditions that are calculated 
from the user input variables Moo, aioo, and /3oo- 


P = PPoo , U = UCoo, V = VC 00 , W = WCoo, p = p / 0 oo C^ o 

G = G Poo Cqq , C = C Coo i T = T Too = C Xqo 
Poo 1) Poo 1 /Tj Gqo 1) Tqo 1 

Uqo — Moo COSOoo COS /3 qo > Voo — M^oo sin/Jooi Woo VIqc 


Moo sin ctoo cos/3oo (5) 


( 2 ) 

(3) 

(4) 


Type Specify Extrapolate Update 


Outflow Moo, [Moo > 0], «oo, Poo None All 


All 



U 


OO 


Moo COSOoo COS/3oo 
-Moo sin/3oo 
Moo sin (x oo cos /Joo 


( 6 ) 
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The right-state vector is: 


Qr = 


Poo 



_ Poo 


(7) 


2.2 Riemann Invariant Boundary Condition (riemann) 

Type Specify Extrapolate Update 

Inflow/outflow Mqo [Moo > 0] Entropy (outflow) p , p 


The Riemann invariants correspond to the incoming R _ and outgoing R + characteristic 
waves. The invariants determine the locally normal velocity component and the speed of 
sound. The entropy, s, and the speed of sound, c, are used to determine the density and 
pressure on the boundary (see Figure 1.) The incoming Riemann invariant uses the far-field 
conditions, q Q = [p x , Uoo,Poo] ■ 


where 


R + = U; + 


2d 

7 - 1 ’ 


R" =U Q 


2c 0 
7 — 1 


Ui = Uj-n, 
U 0 — U 0 -n, 



Ci 


U; = [u i ,v i ,w i ] T , 
U 0 = [u 0 , v Q , w q ] t 


c? = 


7Pi 

Pi 

= 7Po 
Po 


(8) 

(9) 

( 10 ) 

( 11 ) 


If the flow at the boundary is locally supersonic leaving the domain (M, > 1), then no 
incoming characteristic waves exist; thus, R _ is set equal to 


R~ = Ui - 



(12) 


Similarly, if the flow is supersonic entering the domain, then no outgoing characteristic 
waves exist and R + is set equal to 


R+ = U 0 + (13) 

7-1 

A velocity, Ub, and the speed of sound, Cb, at the boundary are the sum and difference of 
the invariants: 


U b = ^ (R + + R”) , c b = 4 ( 7 - 1) (R+ - R-) (14) 

The velocity that is imposed on the boundary depends on the local direction of flow. If 
the sign of U_l is positive, then the flow is exiting the computational domain and the 
entropy is extrapolated from the interior and is used to update the density at the boundary. 
Conversely, a negative Uj_ indicates that the flow is entering the computational domain and 
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the free-stream entropy is used. Summarizing, the velocity and entropy on the boundary 
are calculated from the following equations. 


U b = 


Sb = < 7Pl 


Ui + (U — Uj)n, if U > 0 (outflow) 
U 0 + (U — U 0 )n, if U < 0 (inflow) 

if U_l > 0 (outflow) 
if Uj_ < 0 (inflow) 


The density and pressure on the boundary are then calculated as follows. 


* - j 

\7 s b/ 


1 / 7-1 


Pb = 


PbC 2 h 


The right-side vector is 


(15) 

(16) 


(17) 


Pb 


Q r = 


U b 

.Pb. 


(18) 


2.3 Outflow Mach-Number Boundary Condition (subsonic.outf low_mach) 

Type Specify Extrapolate Update 

Outflow M [0 < M < 1] T p, p 


The Mach number at the boundary, M set , is specified for the entire boundary, and the flow 
is assumed to be adiabatic and isentropic. The local acoustic speed is determined from 
the temperature of the interior state, q ; , with the local speed of sound determined using 
Eq. (19). 

Ci = Vrl, Ti = — (19) 

Pi 

For each point on the boundary, the Mach number, Mj, is used to determine a new total 
pressure, p ti b- 

M; = M Ui=Ui-n (20) 

Ci 

■y 

7-1 

(21) 


Pt,b = Pi 


l+-( 7 — 1)M? 


The static pressure, p, is updated using the new total pressure and the set Mach number, 
M set 


P = Pt,b 


1 + \(.1 - 1)ML 


1 

7-1 


(22) 


If the flow is locally supersonic as a result of some transient flow condition, then the boundary 
condition will reset the local static pressure to the local total pressure. If the specified Mach 
number is supersonic, then an extrapolation boundary condition should be used. Thus, 


Pb 



if Mj < 1 (subsonic) 
if Mi > 1 (supersonic) 


(23) 
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The right-state state is written as: 


Qr = 


7Pb/ T i' 

Ui 

L Pb . 


(24) 


Note: Using this condition is inconsistent if the boundary face adjoins a viscous surface. 
The Mach number is fixed over the entire boundary and precludes the presence of a velocity 
gradient that would normally be present on the viscous no-slip surface. A constant pressure 
boundary condition, such as the subsonic outflow boundary, should be used in this situation. 


2.4 Pressure Outflow Boundary Condition (back_pressure) 

Type Specify Extrapolate Update 

Outflow p/poo [p/Poo > o] U, T p 


For a subsonic outflow boundary condition, the static pressure ratio, p se t/Poo) is specified, 
while the velocities and temperature are extrapolated. The boundary flow is assumed to be 
adiabatic and isentropic. The density is updated from the extrapolated temperature and 
the requested static pressure. If the flow is supersonic at this boundary, then all quantities 
are extrapolated. The pressure is calculated as 


Pb 


Pset 5 

Pi, 


if Mi < 1 (subsonic) 
if Mi > 1 (supersonic) 


and the right-side vector is 


Qr = 


7Pb/T( 

Ui 


Pb J 


T; 


TPi 

Pi 


(25) 


(26) 


2.5 Subsonic Outflow Boundary Condition (subsonic_outf low_p0) 

Type Specify Extrapolate Update 

Outflow p/poo [p/poo > o] U, T p 


The manner in which this boundary specifies the static pressure ratio is the same as that 
presented in section 2.4 but with different implementation details. If any reverse flow (i.e., 
flow into the computational domain) occurs, setting the static pressure at the boundary 
is numerically ill-posed. This boundary condition will explicitly set the flow to exit the 
domain. The flow is also forced to remain subsonic by setting the local static pressure to 
the local total pressure if the local Mach number is greater than one. The velocities, Ui, and 
the temperature, T, = 7 Pi/pi, are extrapolated from the interior solution. The pressure is 

p b = j Pset ’ if Mi < 1 (subsonic) 

[pt, if Mi > 1 (supersonic) 
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and the velocity is 


U b = 


Ui, if U_l > 0 (outflow) 
|Uj|n, if Uj_ < 0 (inflow) 


U_L = Ui -n 


The right-side vector is 


<?r = 


7Pb/ T i' 

Ub 

Pset 


2.6 Mass Flow Out Boundary Condition (massf lux_out) 

Type Specify Extrapolate Update 

Outflow m [m > 0] U , T p, p 


(28) 


(29) 


This boundary condition allows for specification of the mass flow out of the computational 
domain. Adiabatic flow through the boundary is assumed, and the method iteratively 
modifies the static pressure for the entire boundary to obtain the requested mass flow 
condition. The choice of equations solved for this boundary condition is from the CFD code 
VULCAN [10]. 

The boundary values of momentum thrust Eq.(30), pressure force Eq.(31), and mass flow 
Eq.(32), are calculated by using the following integrals, where the domain of the integration 
is the entire boundary face: 


F m = 

f p f U • n'j dA, 

J boundary ^ / 

(momentum thrust) 

(30) 

F P = 

/ prfA, 

J boundary 

(pressure force) 

(31) 

m = 

f p ^U • n"j dA, 

J boundary ^ / 

(mass flow) 

(32) 

The static pressure is updated with a form of the 

continuity equation: 



Pb = 


F m (5m + F r 


(5m =11 — 


ihset \ 


- 


(33) 


The velocity and the temperature are extrapolated from the interior computational domain 
so that the right-state vector is: 


Qr = 


7Pb/ T i 

Ui 

L Pb J 


(34) 


2.7 Subsonic Inflow Boundary Condition (subsonic_inf low_pt) 

Type Specify Extrapolate Update 


Inflow Pt/Poo, Tt/Too, n p, R+, H t 


P, u 


7 


The user can specify the direction of the velocity at the boundary either by the flow angles 
(a se t and /3 se t ) relative to the global coordinate axis, or as normal to the boundary, rfy 


cos a set COS p set i - sin f3 set j + sin a set sin /3 set k 
fib 


(35) 


The right-side vector, q R , is determined from the outward propagating invariant, R + , and 
a statement of the total enthalpy, H t , at an element face on the boundary. The formulation 
of this boundary condition is the same as that used in [10], which assumes that the flow 
through the boundary is adiabatic and isentropic. Following these assumptions, we write 


H ti = » 


if 7 A I 1 / 2 I 2 , 2\ 

■ J + 2 ( U i +^i +W i) 


R+ = -Ui- 


2ci 

7 - 1 


(36) 

(37) 


Because the flow is adiabatic, that is, the total enthalpy is conserved across the boundary, 
we can state that 


Ht = 


b + ^l 


7-12 

Then, by extrapolating R~ to the boundary, we obtain 

2cb 


R + = -U b - 


7 - 1 


By combining Eqs. (39) and (38), we obtain 


Ht = 


7-1 


1 ( R + 


_2cb_ 

7-1 


(38) 


(39) 


(40) 


Equation (40) is solved for Cb, which is the sonic speed at the boundary, by rewriting it as 
a quadratic equation of the form 


1+ (7-l)J 

The solution has the form 


c 2 h + 2R+c b + 7 2 1 (r + 2 - 2H t ) = 0 


b \/fe 2 — 4ac 


Cb± 2 a ± 2a 
where a, 6, and c are the coefficients of the quadratic equation (Eq. 41): 


CL — 1 -|- 


(7-1)’ 


b= 2R+, c = 


7-1 


(r+ 2 - 2H t ) 


(41) 


(42) 


(43) 


The physically consistent result is the larger of the two roots and is chosen to update the 
sonic velocity at the boundary: 


c b = max(c b+ ,c b _) (44) 

The updated inflow velocity and the Mach number at the boundary are then computed by 
using 

u = - R+, Mb = — (45) 

7-1 c b 
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The static pressure and temperature are calculated from the isentropic relations 

7 7 ~ 1 

Pb = Pt,set (l + ^ M b) " _1 , T b = T t , set 7 (46) 

The primitive variables that are determined from the input boundary conditions are placed 
in the right state q R as 


<?R = 


Pb/i?T b 

Un 

Pb 


2.8 Mass Flow In Boundary Condition (massf lux_in) 

Type Specify Extrapolate Update 

Inflow m [m > 0], Tt/Too p p 


(47) 


Mass flow into the computational domain through this boundary condition is updated by 
adjusting the static temperature, T, at the boundary through a form of the energy equation: 


u 2 

CpT H — — CpTt lSe t 


(48) 


The variable, T t , S et, in Eq. (48) is the user-specified total temperature at the boundary, 
and the flow is assumed to be isentropic and adiabatic through the boundary face. An 
expression for the velocity is derived by dividing the specified mass flow by the density and 
can be written as 



P = 


A 


boundary 


p dA 


(49) 


A quadratic equation of static temperature is written by combining Eq. (49) with Eq. (48) 
and using a thermally perfect gas assumption for density in terms of the static pressure and 
temperature, p = p/i?T: 


2 ( pA ) T " CpT c P Tt ’ se ‘-° 

The solution to the resulting equation is a quadratic equation 

rp b \/b 2 - 4ac 

2 a 2a 

where the coefficients of the quadratic equation are 

1 / mi?.\ 2 

a = 2 ( J 1 0 = C P’ C= c pTt,set 

The larger root is taken to update the static temperature at the boundary as 


T b = max(T + , T_) 


(50) 


(51) 


(52) 


(53) 
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To lessen transients during the solution startup, the user can ramp the mass inflow, hi, from 
zero up to the specified amount by using the iteration parameter f low_mf lux_ramp as 


= min 1 


iteration \ 
f lowjnf luxjramp J 


m se t 


(54) 


The continuity equation determines the mean velocity at the boundary, and the density is 
updated from the updated static temperature as 


The right-state vector is: 


Tb hi ram p 

7 P A 


Qr = 


Pi 

-Un 

PiTb/7 


(55) 


(56) 


2.9 Supersonic Inflow Boundary Condition (fixed inflow) 

Type Specify Extrapolate Update 

Inflow p, U, p None All 


Pressure, density, and velocity (U se t) are specified for the supersonic inflow boundary con- 
dition. The velocity vector is forced to be normal to the boundary. The right-state vector 
is: 


<?R = 


Pset 
I U set | n 
Pset 


(57) 


3 Computational Method 

3.1 FUN3D Code 

FUN3D is an unstructured three-dimensional, implicit, Navier-Stokes code. Roe’s flux dif- 
ference splitting [11] is used for the calculation of the explicit terms. Other available flux 
construction methods include HLLC [12], AUFS [13], and LDFSS [14]. The default method 
for calculation of the Jacobians is the flux function of van Leer [15], but the method by 
Roe and the HLLC, AUFS and LDFSS methods are also available. The use of flux limiters 
are mesh and flow dependent. Flux limiting options include MinMocl [16] and methods by 
Barth and Jespersen [17] and Venkatakrishnan [18]. Other details regarding FUN3D can 
be found in Anderson and Bonhaus [6] and Anderson et al. [19], as well as in the extensive 
bibliography that is accessible at the FUN3D Web site, http://fun3d.larc.nasa.gov. 


3.2 Boundary Element Discretization 

Discretization of the computational volume consists of any combination of tetrahedra, hex- 
ahedra, prisms, and pyramids. The faces of the volume elements are either triangles or 
quadrilaterals. Therefore, the boundaries consist only of either triangular or quadrilateral 
faces or combinations of both types of face elements. An example of a boundary face that 
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Figure 3. Perspective view of triangular boundary element in Y-Z plane. 
Table 2. Boundary Element Node Weight Factors 

Topology Wi W 2 W 3 

Tetrahedral 6/8 1/8 1/8 

Hexahedral, prism, pyramid 10 0 


consists of triangular face elements is shown in the perspective sketch in Figure 3. A space 
of subelements is formed to calculate the fluxes at the boundary. Three subelements are 
formed from each of the corner nodes (1, 2, and 3) of the triangular face element. The area 
and the orientation of the quadrilateral subelements are determined by using the center of 
the element , P c ; the left and right midpoints, P 1 and P r ; and the main node, P m . The 

normal vector is the cross product of the two vectors that are formed from the main and 

> > > > 

center nodes, P m P c , and the left and right midpoints, P;P r , so that if = P m P c x P;P r . 

The interior, or left state, for the subelement is calculated from an area weighted average 
of the interior state nodal values, q 1 , q 2 , and q 3 . The values of the area weights, W, change 
with the element topology and are listed in Table 2. A detailed derivation of the boundary 
node weighting for tetrahedral meshes can be found in appendix A. 

<?l = VViQi + W 2 q 2 + W 3 q 3 (58) 

A flux function calculates the contribution to the solution residual for the subelement from 
the interior and right-side (i.e., boundary condition ) states. The inviscid flux is added to 
the solution residual at each solution time step. Roe’s approximate Riemann solver was 
used for the boundary flux contributions to the residual. 


4 Test Case Descriptions 

The three geometries that are described in Sections 4.1 through 4.3 are used to demonstrate 
the inflow and outflow boundary conditions that were described in Section 2. The geometries 
are representative of a typical physical situation that could be modeled with each of the 
different boundary conditions. Static pressure and mass flow outflow subsonic conditions are 
applied to the outflow boundary of the bell-mouth geometry, Section 4.1. Total pressure- 
total temperature conditions and mass flow inflow conditions are applied to the inflow 
boundary of the American Society Of Mechanical Engineers ( ASME ) nozzle geometry, 
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Section 4.2. Fixed primitive variable conditions that simulate a supersonic inflow are applied 
to the inflow boundary of the diffuser geometry, Section 4.3. 


4.1 Generic Bell-Mouth 

The generic bell-mouth configuration is a geometry that transitions an external flow to 
an internal flow over a range of conditions without flow separation [20]. The bell-mouth 
geometry is routinely used at experimental facilities as a flow-conditioning device for testing 
inlets or propulsion simulators. Figure 4(a) shows a perspective view of the geometry with 
the symmetry planes, the external outflow plane, and the external surface of the bell-mouth. 
A detailed view of the external bell-mouth surface from the same perspective view is shown 



Z-symmetry 
plane 


External 


Y- symmetry 
plane 


(a) Perspective view of the geometry. 



m 




Z-symmetry S8| 
plane 


External 


outflow plane 


Bell mouth 


mk 


Flow I 
direction S 


Y-symmetry 

plane 


(b) Detailed view of external bell-mouth ge- 
ometry. 


Figure 4. Bell-mouth geometry with representative surface mesh. 


in Figure 4(b). The direction of the flow is from left to right. The radius of the inlet is 1.0, 
and the far-held boundary is 30 inlet radii upstream, as shown in Figure 5(a). The shape 




(a) Symmetry plane view. 


(b) Geometric definition of the bell-mouth. 


Figure 5. Bellmouth geometry details. 


of the bell-mouth consists of an elliptic contraction section with a 9:1 capture area. (See 
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(a) Perspective view of the nozzle. (b) Geometric definition of the nozzle. 

Figure 6. ASME nozzle geometry. 

Figure 5(b).) The bell fairs to a flat external wall via a circular-arc lip. For this study, 
the geometry is axisymmetric and is modeled with a quarter-plane symmetric unstructured 
tetrahedral mesh. The mathematical equation that describes the elliptic portion of the 
bell-mouth is 

(a - h e ) + (y- K) = he = 0 0) ke = 3.0, a e = 4.5, b e = 2.0 (59) 

and the equation that describes the circular-arc fairing is 

(x - h c ) 2 + (y - k c ) 2 = R^, h c = -a e + R c , k c = k e , R c = 1.0 (60) 

4.2 ASME Flow Calibration Nozzle 

The ASME nozzle is a flow standard geometry that is used for calibration of full-pipe 
flows in experimental facilities [21] and serves as an excellent configuration for validation of 
internal flow calculations. A perspective view of the mesh is shown in Figure 6(a), and a 
streamwise cross-sectional view of the nozzle is shown in Figure 6(b). The geometry consists 
of two constant-area cylindrical sections joined via two elliptically contoured sections. The 
upstream plenum section has a cross-sectional radius of 5 in., and the downstream nozzle 
section is a constant-area pipe with a radius of 2.5 in. The plenum entrance is at station 
x = —10.0, with the exit at station x = 8.0. The fairing tangent point between the two 
elliptical contours is at x ~ 0.54. All of the computations were accomplished by using 
quarter-plane symmetric tetrahedral meshes. The equation that defines the nozzle wall 
contour from the middle tangent point to the downstream tangent point is 

Or — h n ) + (j/-k n ) h n = 5.0, k n = 2.5, a n = 5.0, b n = 2.5 (61) 

a n b n 

and the equation that defines the nozzle wall from the upstream tangent point to the middle 
tangent point is 

(x ~ 2 hp) + ~ 2 kp) = 1, h p = -3.0, k p = 2.5, a p = 3.99, b p = 2.9 (62) 
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Inflow 

plane 


(a) Perspective view of the diffuser. 



(b) Geometric definition of the diffuser. 


Figure 7. Conical diffuser geometry. 
Table 3. Test Case Mesh Densities 


An 

Bell-mouth 

ASME nozzle 

Diffuser 


tetrahedra/nodes 

tetrahedra/nodes 

tetrahedra/nodes 

1.0 

51,078/10,253 

50,309/10,202 

25,484/5,670 

0.8 

99,125/19,326 

97,517/19,129 

52,499/10,992 

0.6 

235,111/44,370 

231,510/43,760 

132,693/26,452 

0.4 

775,258/141,372 

764,440/139,535 

433,863/82,184 

0.2 

6,648,653/1,168,680 

6,011,766/1,059,652 

2,271,735/491,521 


4.3 Supersonic Diffuser 

The third test geometry is a conical diffuser with a wall divergence angle of 3 degrees. 
The increasing cross-sectional area of the duct maintains the supersonic flow condition 
downstream of the supersonic fixed inflow boundary. A perspective view of the geometry 
is shown in Figure 7(a). A streamwise cross-sectional view of the geometry is shown in 
Figure 7(b). 


5 Discussion 

The unstructured meshes with tetrahedral cells were generated using the software package 
VGRID [22]. The number of nodes for the cases is changed through the use of the global 
parameter, ifact. The ifact parameter controls the density and number of cells, N, by 
increasing or decreasing the strength of the sources within the mesh generation data file. 
The total number of cells (and nodes) for each case are listed in Table 3. 

An equivalent mesh size can be related to the mesh density (or number of degrees of 
freedom) N. The equivalent mesh size should decrease with an increase in the number of 
cells. In three dimensions, the equivalent mesh size, h]\j, should tend to zero as N -1 / 3 . The 
equivalent mesh size can also be related to a characteristic distance defined in terms of the 
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0 Bell mouth 

— — Q — ASME nozzle 
€> — Diffuser 



Equivalent mesh size, h N 


Figure 8. Mesh characteristics. 
Table 4. Test Cases 


Boundary condition 

Boundary name 

Geometry 

Specified conditions 

Pressure outflow 

backjpressure 

Bell-mouth 

p/Poo 

Pressure outflow 

subsonic_outf low_pO 

Bell-mouth 

p/Poo 

Mach number outflow 

subsonic_outf low_mach 

Bell-mouth 

M 

Mass flow outflow 

massf lux_out 

Bell-mouth 

m 

Pressure subsonic inflow 

subsonic_inf low_pt 

ASME nozzle 

Pt/Poo, Tt/Too , a, (3 

Mass flow inflow 

massf lux_in 

ASME nozzle 

hi, Tt/Too 

Supersonic inflow 

f ixed_inf low 

Diffuser 

P, U, p 


local control volume size, h v = ||V 1</3 ||, where the norm used in this study is the L,*, norm. 


A statement of the asymptotic order of error can only be made when consistently refined 
meshes are used. Consistent mesh refinement is purely geometric occurring when a linear 
relationship exists between the equivalent mesh size based on degrees of freedom ,1 in, and 
the equivalent mesh size based on the characteristic distance, hy. See Thomas et al. [23] 
for additional discussion on consistent refinement. The similarity of the mesh families for 
each configuration is shown in Figure 8. A consistently refined mesh is shown as a dashed 
line in the Figure. All three geometries display a linear trend between the equivalent mesh 
measures, h^ and hy. 

The free-stream Mach number, unless otherwise noted, was 0.2 and the Courant-Friedrichs- 
Lewy (CFL) number, unless otherwise noted, was ramped from 1 to 100 over 250 iterations. 
Flux construction was performed by using the approximate Riemann solver by Roe and 
without any flux limiting. The boundary condition test cases are listed in Table 4. The 
bell-mouth was used to test the static pressure and the mass flow out boundary conditions. 
The ASME nozzle was used to test boundary conditions with the total temperature and ei- 
ther total pressure or mass flow set at an inflow boundary. The diffuser geometry was used 
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Figure 9. Mach number contours for the bell-mouth geometry. 


to test the supersonic fixed inflow boundary condition. In the following sections, solution 
convergence histories for each of the mesh densities are shown for each boundary condi- 
tion. Also, the convergence of the error in the boundary value is shown for each boundary 
condition. The error in the boundary condition is the difference between the user-defined 
boundary condition and the area-weighted average of the boundary solution condition from 
the code. The area-weighted average for a variable, <j>, is designated by an overbar <f> and is 
calculated as / boundary / boundary dA. 

The area- weighted average values for p, p, p t , and mass flow are calculated by the 
summation equations, Eqs. (63). The average velocity, speed of sound, and Mach number 
are calculated by Eqs. (64). 

p = \ p 5A ’ p = x p<5A ’ Ft = x (p+^u 2 V A ( 63 ) 

boundary boundary boundary ^ / 

m = pU_i_<5A, U = m/pA, c = yj yp /p, M = U/c (64) 

boundary 


where JA is the area of an individual element on the boundary face and A is the boundary 
area, A = ^boundary ^ A - Th e error in the boundary condition for the parameter <f> is then 


Error ^ = (j) set - (j) 


(65) 


5.1 Bell-Mouth Calculations 

The bell-mouth geometry was used to test the fixed subsonic Mach number, fixed static 
pressure, and mass flow outflow boundary conditions. Surface Mach number contours on 
and around the bell-mouth for an outflow mass flow condition of 0.25 and a free-stream 
Mach number of 0.2 are shown in Figure 9. For this set of conditions, the flow stagnation 
point is slightly below the leading edge of the bell-mouth, as evidenced by the region of zero 
Mach number contours. The flow accelerates externally around the bell mouth leading edge 
to exit the computational domain further downstream. 
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5.1.1 Static Pressure Outflow (back_pressure) 

Type Specify 

Outflow p/poo = 0.90, Moo = 0.2 




(a) Residual history. (b) Boundary value error history. 



10 1 ‘ ‘ — 

0.1 0.3 0.5 0.7 0.91.1 


Equivalent mesh size, h N 
(c) Boundary value error property. 


Figure 10. Static pressure ratio boundary condition, back_pressure. 

The solution residual history for the fixed static pressure (also known as a back pressure) test 
is shown in Figure 10(a). Iterative convergence was achieved for all mesh levels for the set 
condition of p/poo = 0.9. The boundary value error typically reached the final level for each 
equivalent mesh size in fewer than 400 iterations, Figure. 10(b). The final boundary value 
error is typically achieved before interative convergence of the solution residual. The change 
in boundary value error with equivalent mesh size after iterative convergence is shown in 
Figure 10(c). The error in the set static pressure converges as the square of the equivalent 
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mesh size. The second-order convergence property in the error is indicated by the slope of 
the dotted line. 

5.1.2 Subsonic Static Pressure Outflow (subsonic_outf low_pO) 

Type Specify 

Outflow p/poo = 0.90, Mqo = 0.2 




(a) Residual history. (b) Boundary value error history. 



10 1 ‘ ‘ 

0.1 0.3 0.5 0.7 0.91.1 


Equivalent mesh size, h N 
(c) Boundary value error property. 


Figure 11. Subsonic outflow boundary condition, subsonic.outf low_p0. 

The residual history for subsonic outflow static pressure calculations is shown in Fig- 
ure 11(a). Similar to the back pressure boundary, all mesh levels achieve iterative con- 
vergence for the set condition of p/poo =0.9. The boundary value error achieved the final 
value for each mesh level typically before 400 iterations, Figure 11(b). The change in bound- 
ary value error with equivalent mesh size after iterative convergence is shown in Figure 11(c). 
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The error in the set static pressure converges as the square of the equivalent mesh size. The 
second-order convergence property in the error is indicated by the slope of the dotted line. 


5.1.3 Mach-number outflow (subsonic_outf low_mach.) 

Type Specify 

Outflow M set = 0.90, Moo= 0.2 



(a) Residual history. (b) Boundary value error history. 



10 1 ‘ ‘ 

0.1 0.3 0.5 0.7 0.91.1 


Equivalent mesh size, h N 
(c) Boundary value error property. 


Figure 12. Mach-number boundary condition, subsonic_outf lowjnach. 

The Mach-number outflow boundary condition adjusts the static pressure to obtain the 
specified Mach number. The residual history for the fixed subsonic Mach-number outflow 
boundary condition test is shown in Figure 12(a). Iterative convergence was achieved for 
all mesh levels for the set condition of M se t = 0.9. As shown in Figure 12(b), the error in 
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the Mach number reaches the final value within several hundred iterations for this test case. 
The error in the set Mach number converges as the square of the equivalent mesh size. The 
second-order convergence property in the error is indicated by the slope of the dotted line 


(Figure 12(c)). For this configuration, M se t= 
static pressure ratio of p/poo = 0.61. 

5.1.4 Mass flow outflow (massf lux_out) 

Type 

Outflow m = 0 



(a) Residual history. 


0.9 is approximately equivalent to a inflow 


Specify 

.25, Moo = 0.2 



(b) Boundary value error history. 



0.1 0.3 0.5 0.7 0.91.1 

Equivalent mesh size, h N 

(c) Boundary value error property. 


Figure 13. Mass flow out boundary condition, massf lux_out. 

The residual history and the boundary value error history for the fixed mass flow outflow 
calculation are shown in Figure 13(a). Iterative convergence was achieved for all mesh levels 
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for the set condition of m = 0.25. As shown in Figure 13(b), the error in the mass flow 
reaches the final value within several hundred iterations for most cases. The relaxation of 
the static pressure to set the mass flow took longer for the equivalent size mesh Iin = 0.2 
as a result of a lower damping of the pressure oscillations in the flow field. The error in 
the set mass flow converges as the square of the equivalent mesh size. The second-order 
convergence property in the error is indicated by the slope of the dotted line (Figure 13(c)). 
The average values for the Mach number and the static pressure ratio at the boundary for 
this mass flow setting are approximately 0.33 and 0.95, respectively. 


5.2 ASME Nozzle Calculations 


Two inflow boundary conditions are evaluated in this section using an ASME flow calibration 
nozzle. Figure 14 shows Mach-number contours for choked-flow conditions in the nozzle. The 
direction of flow, from left to right, is denoted by the symbol U and the vector arrow. Typical 
combinations of inflow and outflow boundaries for an internal flow configuration includes 
(a) total pressure-total temperature inflow with static pressure outflow; (b) total pressure- 
total temperature inflow with mass flow outflow; or (c) mass flow inflow and extrapolation 
outflow. Choice (c) would work for choked flows where the exit Mach number is equal to 
or greater than 1. The total pressure-total temperature inflow boundary condition is the 
most appropriate inflow condition to use when the performance of the nozzle (e.g. discharge 
coefficient, thrust ratio, or thrust coefficient) is to be determined. This set of conditions does 
not presuppose the mass flow of the nozzle because the performance is a result of geometric 
characteristics. 



Figure 14. Mach-number contours for the ASME nozzle test case. 


5.2.1 Pressure Subsonic Inflow (subsonic.inf low_pt) 


Type Specify 

Inflow pt/poo = 1-6, Tt/Too = 1.0, a = (3 = 0.0 
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(a) Residual history. (b) Boundary value error history. 



0.1 0.3 0.5 0.7 0.91.1 

Equivalent mesh size, h N 

(c) Boundary value error property. 

Figure 15. Subsonic inflow (p t ) boundary condition, subsonic_inf lowjpt. 


Iterative convergence, shown by the residual history plotted in Figure 15(a), was achieved 
for all mesh levels for the set condition of pt/poo = 1.6. As shown in Figure 15(b), the error 
in the total pressure reaches the final value within several hundred iterations for all grids. 
The error in the set total pressure converges as the square of the equivalent mesh size. The 
second-order convergence property in the error is indicated by the slope of the dotted line 
(Figure 15(c)). The average values for the Mach number and the mass flow ratio at the 
boundary for this total pressure/total temperature setting are approximately 0.14 and 4.4, 
respectively. 
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5.2.2 


Mass Flow Inflow (massf lux_in) 

Type Specify 


Inflow rii = 1.0, Tt/Too = 1.0 (unchoked) 
Inflow rii = 10.0, Tt/Too = 1.0 (choked) 




(a) Residual history. (b) Boundary value error history. 



10 1 1 1 

0.1 0.3 0.5 0.7 0.91.1 


Equivalent mesh size, h N 
(c) Boundary value error property. 


Figure 16. Unchoked mass flow boundary condition, massflux_in. 

The residual history and the boundary value error history for the mass flow inflow cases 
are shown in Figure 16 for the unchoked nozzle flow and Figure 17 for the choked nozzle 
flow condition. Iterative convergence was achieved for all mesh levels for both set conditions 
of hi = 1.0 and rii = 10.0, as shown in Figures 16(a) and 17(a) respectively. The error 
in the mass flow reaches the final value in about 1000 iterations for all cases, as shown 
in Figures 16(b) and 17(b). The boundary value error converges approximately as the 
square of the equivalent mesh size. The second-order convergence property in the error is 
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Residual (p) 




(a) Residual history. (b) Boundary value error history. 



0.1 0.3 0.5 0.7 0.91.1 

Equivalent mesh size, h N 

(c) Boundary value error property. 


Figure 17. Choked mass flow boundary condition, massflux_in. 
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indicated by the slope of the dotted line (Figures 16(c) and 17(c)). The average values for 
the total pressure ratio at the inflow boundary for the unchoked and choked nozzle mass 
flow conditions are approximately 1.03 and 3.5, respectively. 

5.3 Supersonic Diffuser Calculations 

A solution using the supersonic fixed inflow boundary condition at the inflow plane of the 
diffuser is shown in Figure 18. An extrapolation boundary condition was used at the outflow 
plane. The flow is from left to right, as indicated by the arrow. The flow accelerates from 
approximately Mach 1.3 at the inflow to Mach 2.3 at the outflow. 



plane 


Figure 18. Mach-number contours for the diffuser test case. 


Type Specify 

Inflow p= 1.2, U = (1.0, 0., 0.), p = 1.0/1. 4 


The residual history and the boundary value error history for the fixed supersonic inflow 
case are shown in Figure 19. Iterative convergence was achieved for all mesh levels for the 
set conditions within a few hundred iterations. As shown in Figure 19(b), the error in the 
set velocity achieves the final value for each mesh in less than 100 iterations. The error in 
the set velocity converges as the square of the equivalent mesh size for the three smallest 
equivalent mesh sizes. The second-order convergence property in the error is indicated by 
the slope of the dotted line (Figure 19(c)). 

6 Summary 

Boundary conditions that allow subsonic and supersonic flow into and out of the computa- 
tional domain have been derived and implemented in the computational method FUN3D. 
The boundary conditions are enforced through determination of the flux contribution to 
the solution residual. Roe’s approximate Riemann solver is used to construct the fluxes at 
the boundary faces. These boundary conditions were implemented in an implicit manner in 
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(a) Residual history. (b) Boundary value error history. 



0.1 0.3 0.5 0.7 0.91.1 

Equivalent mesh size, h N 

(c) Boundary value error property. 


Figure 19. Supersonic inflow boundary condition, f ixed_inf low. 
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the FUN3D CFD code and were verified for three generic test geometries. Iterative residual 
convergence was achieved at all mesh-density levels for all of the conditions that were used in 
this study. The user-requested condition was achieved in all cases, and error in the bound- 
ary condition value decreased in a second-order manner with an increase in mesh density. 
Boundary value error convergence occurs before iterative convergence and depending on the 
application, boundary values were within 0.5 percent of the specified parameter value with 
only a residual reduction of 5 orders of magnitude. 
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7 Nomenclature 

Roman letters 
A area 

c speed of sound 

c p specific heat at constant pressure 
H enthalpy 

Iin equivalent mesh size based on degrees of freedom, N -1 / 3 

h v Li norm of V -1 / 3 , |Li| = ||V _1 / 3 ||i 

i,j , k Cartesian unit vectors axis in physical space 

m mass flow 

M Mach number 

n normal vector 

n unit normal vector, n/|n| 

N number of nodes in mesh 

p static pressure 

R real gas constant 

R Riemann invariant 

s entropy 

T temperature 

U average velocity 

U velocity vector, [m,u,u>] T 

U velocity magnitude 

u,v,w Cartesian velocity components 
V primal volume of elements in mesh 

x, y, z Cartesian coordinates in physical space 

Subscripts 

b boundary state 

00 free-stream state 

1 internal state 

L left state 
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o 


outer state 
_L normal to boundary element 

R right state 

set user-requested condition 

t total conditions 

Conventions 

!F boundary element flux vector 

q primitive variables state vector 

1Z solution residual vector 

Symbols 

a u-w angle of velocity 

f3 u-v angle of velocity 

7 ratio of specific heats, 7 = 1.4 (air) 

A eigenvalue 

p density 
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Appendix A 


Node Weighting for a Boundary Tetrahedral 

The tetrahedral element that is shown in Figure A 1 has four nodes with an associated 
area-weighted normal, hi, 112, 113, and 114, opposite each node. The average value of the 



Figure Al. Tetrahedral element. 


variable qi at the face centers is one-third the sum of the three 
assume that a constant gradient exists across the element. 

corners of that face if we 

Ql = ^ (*?node2 Qnode3 H - Qnode^) 

(Al) 

Q 2 = 7^ ((/nodel H" Qnode3 H - Qnode^) 

(A 2 ) 

Q 3 = ^ (^nodel Qnode2 H - Qnode^) 

(A 3 ) 

Qd = ^ (^nodel H - O'node2 H - Qnode 3 ) 

(A 4 ) 


The Green-Gauss (GG) integral over the tetrahedral element is equal to one-third the sum- 
mation of the product of the average value of qt with each directed face normal, n, . 

/•GG 1 4 

® VqdV = WgA (A 5 ) 

•/tetrahedron *-* i = \ 

The dual of the tetrahedral element shown in Figure A 2 has six faces, with the associated 
normals ne, Sh, and Sr on the boundaries and 1112, ni3, and 1114 facing the tetrahedral face 
opposite node 1 . The finite-volume (FV) integral around the dual volume is written as a 
summation of q times the directed normal of the face. At this point, the influence of the 
nodes on the boundary faces L, R, B, a, 6, and c are unknown. The goal is to make the FV 
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Node 3 


Node 2 


Figure A2. Tetrahedral dual element. 


integral exact when the q n are linear. 


/• FV 

f 

./dual 


VqdV 


9i + 92 _ , 9i + 93 . 9i + 94 _ 

— 2 — n i2 + — 2 — ni3 + — 2 — 1114 

+ (agi + bq 3 + cg 4 ) n L 

+ (091 + bq 2 + cq 4 ) Sr 

+ (091 + bq 2 + cq 3 ) n B 


(A6) 


The FV integral around the dual volume can be equated to one-fourth the GG integral 
around the tetrahedron. 


1 

4 


GG 


S7qdV 


tetrahedron 



(A7) 


The following shows that the left-hand side of Eq. (A7) can be expressed in a form that 
matches the right-hand side of Eq. (A6). Equation (A5) is substituted into the right-hand 
side of Eq. (A7) and expanded as 


1 

4 


r GG 


S7qdV 


/ tetrahedron 


1 1 
43 
1 

12 




(< 7 ini + q 2 d 2 + q 3 n 3 + q 4 n 4 ) 


(AS) 

(A9) 


Add and subtract 2qin 4 from Eq. (A9), and use the equations that relate n 2 , n 3 , and n 4 to 
n 4 and the normals of the dual volume, h 4 2 , n 4 3 , and n 44 (Eq. A10). 


hj = 12n 4i + n 4 


(A10) 


1 

4 


(■GG 


VqdV 


/ tetrahedron 


l 

12 


39ini 


29 in 4 + q 2 (12n 12 + n 4 ) 


+ q 3 (12n 43 + n 4 ) + 94 (12n 14 + ni) 


(All) 
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Split the terms with the hii normal vector as 


1 f GG 1 

- / S7qdV = — 

^ ./tetrahedron 


— 2(7x11! + 6 < 7 2 hi 2 + 6(731113 + 6(741114 
+ 3gini + q 2 ni + g 3 hi + q A n 1 
+ 6g 2 ni2 + 6q 3 hi 3 + 6g 4 hi 4 


(A12) 


Substituting into Eq. (A12) the additional geometric relations that equate the face normals 
of hi, ni 2 , n 13 , and hi 4 with n^, h#, and ng (Eqs. (A13)-(A16): 


ni = -3 (ni 2 + hi 3 + hi 4 ) = 3 (n L + n R + n B ) 
5; 12 = ~ (2h_L +n fi + n B ) 

5 13 = -^(h L + 2h fi + h B ) 

514 = (h L + hR + 2 h s ) 


(A13) 

(A14) 

(A15) 

(A16) 


fGG 


' tetrahedron 


Y 7 it r qi + <72 - . qi + 93 . 91 + 94 _ , 1 (o i 

vqdV = — n i2 H n i3 H n i4 + — (3gini) 


+ 92 
+ 93 
+ 94 


(2 n L + Hr + n B ) + * (n L + n K + n B ) 

- ^ (hi + 2n_R + n B ) + ^ (hi + h# + n B ) 

- ^ (hi + n R + 2he) + i (hi, + h# + n B ) 


Finally, gather terms of similar face vectors: 
1 fGG 


-1 Y 7 7 T / 9i + 92 ^ , 9i + 93 ^ . 9i + 94 _ 

- / VqdV = — n i2 H n i3 H n i4 

^ J tetrahedron AAA 

+ (^ 1 + ^ 3 + H lTi 

+ ( l qi + l q2 + l qi ) SR 
+ [ l qi + h 2 + h 3 ) * B 


(A17) 


(A18) 


A direct comparison of the terms in Eq. (A18) with those in Eq. (A6) shows that a = 3/4, 
b = 1/8, and c = 1/8. 


The author would like to thank Nishikawa Hiroaki for this derivation. 
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